#!/usr/bin/python
# GM_logN.py
'''Calculates stars per square degree.'''

import math
import numpy

def logN(simpars, phypars, appearance):
    for i in simpars:
        cmd = "%s = simpars['%s']" % (i,i)
        exec cmd
    for i in phypars:
        cmd = "%s = phypars['%s']" % (i,i)
        exec cmd
    # Relevant distances
    # 1d
    d2dd = d**2*deltad
    cosb = numpy.cos(ba)
    # 2d
    vol = numpy.multiply( d2dd[:,numpy.newaxis], deltala[numpy.newaxis,:] )
    vol = vol*(math.pi**2*180.**-2) # radians**2 -> degrees**2
    # 0d
    if pencilBeam == False:
        lrange = math.pi
        brange = math.pi/2.
        if (lrange*1.01 < deltala.sum()) or (lrange*0.99 > deltala.sum()):
            print 'Hold your horses! Something is wrong with l in GM_logN.'
        if (brange*1.01 < deltaba.sum()) or (brange*0.99 > deltaba.sum()):
            print 'Hold your horses! Something is wrong with b in GM_logN.'
    if pencilBeam == True:
        lrange = float(deltala[0])
        brange = float(deltaba[0])
    bsum = (deltaba*cosb).sum()
    # Integrate
    N = []
    allskyN = []
    for i in range(0, len(appearance)):
        Ni = numpy.multiply( appearance[i], vol[:,:,numpy.newaxis] )
        Ni = Ni.sum(axis=0) # Straight sum over d
        Ni = Ni.sum(axis=0)/lrange # Weighted average over l
        N.append(numpy.log10(Ni))
        Na = Ni*deltaba*cosb
        Na = Na.sum(axis=0)/bsum # Weighted average over b
        allskyN.append(numpy.log10(Na))
    N = numpy.array(N)
    allskyN = numpy.array(allskyN)
    return (N, allskyN)
